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ABSTRACT 

A  de-shadowing  technique  is  presented  for  multispectral  and  hyperspectral  imagery  over  land  acquired  by 
satellite  /  airborne  sensors.  The  method  requires  a  channel  in  the  visible  and  at  least  one  spectral  band  in 
the  near-infrared  (0.8-1  pm)  region,  but  performs  much  better  if  bands  in  the  short-wave  infrared  region 
(around  1.6  and  2.2  pm)  are  available  as  well.  Five  major  processing  steps  are  employed  :  (1)  the 
calculation  of  the  covariance  matrix  and  zero-reflectance  matched  filter  vector  based  on  atmospherically 
corrected  surface  reflectance  data,  (2)  the  derivation  of  the  unsealed  and  scaled  shadow  function,  (3)  a 
histogram  thresholding  of  the  unsealed  shadow  function  to  define  the  core  shadow  areas,  (4)  a  region 
growing  to  include  the  surroundings  of  the  core  shadow  areas  for  a  smooth  shadow/clear  transition,  and 
(5)  the  de-shadowing  of  the  pixels  in  the  final  shadow  mask  using  the  scaled  shadow  function.  A  radiative 
transfer  code  is  employed  in  the  last  step  to  compute  the  de-shadowed  surface  reflectance  imagery. 
Example  images  of  hyperspectral  airborne  and  multispectral  spaceborne  are  presented. 


1.  INTRODUCTION 


Remotely  sensed  optical  imagery  of  the  Earth's  surface  is  often  contaminated  with  cloud  and  cloud 
shadow  areas.  Surface  information  under  cloud  covered  regions  cannot  be  retrieved  with  optical  sensors, 
because  the  signal  contains  no  ground  component.  In  shadow  areas,  however,  the  ground-reflected  solar 
radiance  is  always  a  small  non-zero  signal,  because  the  total  radiation  signal  at  the  sensor  contains  a  direct 
(beam)  and  a  diffuse  (reflected  skylight)  component.  Even  if  the  direct  solar  beam  is  completely  blocked 
in  shadow  regions,  the  reflected  diffuse  flux  will  remain.  Therefore,  an  estimate  of  the  fraction  of  direct 
solar  irradiance  for  a  fully  or  partially  shadowed  pixel  can  be  the  basis  of  a  compensation  process  called 
de-shadowing  or  shadow  removal. 

Several  de-shadowing  techniques  have  been  described  in  the  literature.  One  approach  employs  geometric 
considerations  to  project  the  cloud  structure  on  the  ground  based  on  the  direction  of  the  incident  solar 
radiation  and  the  estimated  cloud  height  (Simpson  and  Stitt  1998).  Another  approach  uses  the  matched 
filter  concept  (Adler-Golden  et  al.  2002)  where  the  matched  filter  vector  is  calculated  in  terms  of  the 
covariance  matrix  of  the  atmospherically  corrected  reflectance  data. 

This  contribution  also  employs  the  matched  filter  concept  as  one  step  of  the  de-shadowing  method.  The 
application  of  the  matched  filter  vector  to  the  atmospherically  corrected  reflectance  imagery  yields 
unsealed,  positive  and  negative  shadow  abundance  values.  However,  the  re-scaling  of  the  unsealed 
function  into  the  physical  interval  (0,1),  containing  the  fraction  of  direct  illumination,  is  performed  with  a 
new  strategy,  as  detailed  in  the  next  section. 


Richter,  R.  (2005)  De-Shadowing  of  Satellite/ Airborne  Multispectral  and  Hyperspectral  Imagery.  In  Emerging  EO  Phenomenology 
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ORGANIZATION 


All  mentioned  spectral  methods  imply  the  linear  independence  of  pure  spectra.  This  may  not  be  true, 
especially  in  complex  scenes.  So  mis-classifications  do  occur.  As  an  example,  a  fully  illuminated  pixel  of 
a  dark  material  can  be  classified  as  a  shadowed  pixel  of  another  material  or  combination  of  materials.  Still, 
useful  de-shadowed  images  have  been  processed  for  a  variety  of  different  climatic  and  landscape 
conditions. 

2.  DE-SHADOWING  METHOD 

The  method  consists  of  a  sequence  of  processing  steps,  see  Fig.  1.  It  starts  with  a  calculation  of  the  surface 
reflectance  image  cube  pi  —  p( )  ,  where  three  spectral  bands  around  Aj  =  0.85,  1.6,  and  2.2  pm  are 

selected.  These  bands  from  the  near  and  shortwave  infrared  region  are  very  sensitive  to  cloud  shadow 
effects,  because  the  direct  part  of  the  solar  radiation  flux  at  the  ground  level  is  typically  80%  or  more  of 
the  total  downwelling  flux,  see  Fig.  2.  Channels  in  the  blue-to-red  region  (0.4  -  0.7  pm)  are  not  used  for 
the  detection  of  shadow  regions  because  they  receive  a  much  larger  diffuse  radiation  component,  making 
them  less  sensitive  to  partial  shadow  effects.  In  this  contribution  the  ATCOR  code  is  used  for  the 
atmospheric  correction  (Richter  1996,  1998,  Richter  and  Schlapfer  2002,  http://www.rese.ch  ),  which 
employs  a  database  of  radiative  transfer  calculations  based  upon  MODTRAN  (Berk  et  al.,  2000). 


Fig.  1:  Processing  steps  of  de-shadowing  method. 


The  surface  reflectance  is  first  computed  with  the  assumption  of  full  solar  illumination,  i.e.,  the  global  flux 
on  the  ground  consists  of  the  direct  (Edir )  and  diffuse  (Edif)  component.  If  DN  denotes  the  digital  number 
of  a  pixel,  LP  the  path  radiance,  and  x  the  atmospheric  transmittance  (ground-to-sensor)  the  surface 
reflectance  can  be  obtained  as  (Richter  1996)  : 


Pi(x.y) 


(1) 


Here,  d  is  the  Earth-Sun  distance  at  the  image  acqusition  time  in  astronomical  units,  c0  and  Ci  are  the 
radiometric  calibration  coefficients  (offset  and  slope)  to  convert  the  digital  number  into  the  corresponding 
at-sensor  radiance  L,  i.e.,  L  =  c0+  Ci  DN  ,  and  i  is  the  channel  index. 
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The  next  step  is  the  masking  of  water  bodies  and  cloud  areas  with  simple  spectral  criteria  as  detailed  in  the 
following  discussion  section.  Water  pixels  have  to  be  excluded  as  far  as  possible  to  avoid  their  assignment 
as  shadow  pixels.  In  addition,  cloud  pixels  should  not  be  processed.  This  is  achieved  with  the  following 
thresholds: 


p(0.48  pm)  >  30  %  and  p(l .6  pm)  >  30  %  (cloud)  (la) 

p(0.85pm)<  5%  and  p(1.6pm)<  1%  (water)  (lb) 

Step  3  calculates  the  covariance  matrix  C( p' )  where  pf  is  the  reflectance  vector  comprising  only  the 
non-water  and  non-cloud  pixels.  For  each  pixel,  this  vector  holds  the  reflectance  values  in  the  3  selected 
channels  (around  0.85,  1.6,  2.2  pm).  The  matched  filter  is  a  vector  tuned  to  a  certain  target  reflectance 
spectrum  p't  to  be  detected  (Adler-Golden  et  al.  2002) : 


Vmf 


c~V  p;  -  p  v 

rp;-p'/c-Yp;-p7 


(2) 


Here,  pf  is  the  scene-average  spectrum,  without  the  water/cloud  pixels.  Selecting  p't=  0  for  a  shadow 
target  yields  a  special  simplified  form  of  the  matched  filter,  where  the  ?sh?  index  symbolizes  shadow  : 


Vsh=- 


(3) 


The  shadow  matched  filter  vector  is  then  applied  to  the  non-water/non-cloud  part  of  the  scene  and  yields 
the  still  un-normalized  values  O  that  are  a  relative  measure  of  the  fractional  direct  illumination,  also  called 
unsealed  shadow  function  here  : 

0(x,y)  =  V^h{p'(x,y)-p')  (4) 

The  matched  filter  calculates  a  minimum  RMS  shadow  target  abundance  for  the  entire  (non-water/non- 
cloud)  scene.  Therefore,  the  values  of  O  are  positive  and  negative  numbers. 

The  arbitrary,  image-depending  range  of  O  has  to  be  rescaled  to  the  physical  range  from  0  to  1,  where  0 
indicates  no  direct  illumination  (full  shadow),  and  1  means  full  direct  illumination.  The  histogram  of  O  is 
used  to  rescale  the  data.  For  this  purpose,  we  assume  that  the  location  of  the  main  peak  at  Omax  represents 
fully  illuminated  areas.  It  is  based  on  the  statistical  reasoning  that  full  direct  solar  illumination  is  already 
obtained  for  surfaces  of  average  scene  brightness.  Then  the  values  O  are  linearly  mapped  from  the 
unsealed  (Omin,  Omax)  interval  onto  the  physically  scaled  (0,1)  interval,  where  the  scaled  shadow  function 

is  named  Q*  : 


0*  =  Q-^min 


0  -O 

^ max  ^ min 


if  0  <  0 

11  ^  ^  Tt 


0  =  1 


if  0>0n 


(5a) 

(5b) 


The  smallest  value  of  the  scaled  shadow  function  is  0*m-m  =  0,  which  means  no  direct  illumination. 
However,  to  avoid  overcorrection  and  to  cope  with  scenes  containing  merely  partial  shadow  areas,  it  is 
advisable  to  set  0*nin  at  a  small  positive  value  (typically  between  0.05  and  0.10). 
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ORGANIZATION 


In  principle,  the  de-shadowing  could  now  be  performed  with  the  physically  scaled  function  0* ,  which 

represents  the  fraction  of  the  direct  illumination  for  each  pixel  in  the  p f  vector,  i.e.,  the  complete  scene 
without  cloud  and  water  pixels.  However,  since  the  matched  filter  is  not  a  perfect  shadow  transformation, 
it  is  much  better  to  restrict  its  application  to  the  potential,  most-likely  shadow  areas.  This  is  an  important 
processing  step  to  reduce  the  number  of  mis-classifications  or  false-alarms.  If  omitted  it  will  cause  strange 
'shadow'  pixels  scattered  all  over  the  image,  compare  Fig.  2. 


Fig.  2: 

Left:  surface  reflectance  image  of  HyMap  at  Chinchon,  Spain.  Colour  coding:  RGB=878,  646,  462  nm 
channels. 

Center:  standard  shadow  map. 

Right:  map  derived  from  core  shadow  regions. 


Therefore,  the  proposed  method  tries  to  find  the  core  shadow  areas  in  a  scene,  and  subsequently  expands 
the  core  regions  to  obtain  the  final  mask  that  includes  a  smooth  shadow/clear  transition.  The  physically 

scaled  shadow  function  0*  is  then  applied  only  to  the  pixels  in  the  final  mask. 


The  histogram  of  the  unsealed  shadow  function  O  can  be  employed  to  separate  regions  of  low  values  of  O 
from  the  moderate-to-high  values  (see  Richter  and  Muller  2005  for  details).  A  threshold  Ot  can  be  set  in 
the  vicinity  of  the  local  histogram  minimum  and  the  core  shadow  mask  is  defined  by  those  pixels  with 
0(x,y)  <  0T  • 
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Once  the  core  shadow  mask  has  been  defined,  it  is  expanded  to  include  the  surrounding  shadow/clear 
transition  zone  of  100  m  width.  De-shadowing  with  the  scaled  shadow  function  is  then  exclusively 
applied  to  the  pixels  in  this  final  mask.  This  means  the  direct  solar  flux  (Edir  term  in  eq.  1)  has  to  be 
multiplied  with  0*  (x,  y)  : 


Pi(*>y) 


2  {co(i)  +  ci(i)DNi(x>y)}-LPji) 

Ti  [Edirj  0*(x>y)+  Edif,i) 


(6) 


3.  EXAMPLES  OF  DE-SHADOWING  IMAGERY 


Three  selected  examples  show  the  performance  of  the  method  for  a  hyperspectral  airborne  scene  and  two 
multispectral  satellite  scenes. 

The  first  example  of  processing  (Fig.  2,  left)  demonstrates  a  successful  shadow  removal  for  a  situation 
where  all  clouds  are  outside  the  airborne  scene,  so  a  cloud  shadow  algorithm  relying  on  cloud  patterns 
would  fail.  The  imagery  was  recorded  by  the  HyMap  sensor  (Cocks  et  al.  1998)  covering  the  0.45  -  2.45 
pm  spectral  region  with  126  bands. 

Fig.  2  (left)  shows  the  atmospherically  corrected  surface  reflectance  image  prior  to  de-shadowing.  The 
center  of  the  figure  presents  the  standard  shadow  map  employing  all  image  pixels.  The  right  part  is  the 
corresponding  improved  shadow  map  derived  from  the  core  shadow  areas.  Obviously,  the  standard 
shadow  map  contains  a  lot  of  misclassified  pixels,  such  as  roads,  borders  between  different  fields,  and 
others.  The  improved  shadow  map  also  contains  some  small  artifact  areas,  but  it  has  a  much  lower  number 
of  misclassified  pixels.  It  nicely  captures  the  significant  shadow  areas.  Fig.  3  shows  a  zoomed  view  from 
the  top  of  this  scene.  It  presents  reflectance  data  before  de-shadowing,  and  results  after  de-shadowing  with 
the  core-shadow  method. 


Fig.  3  HyMap  scene,  12  July  2003,  Chinchon,  Spain,  RGB  =  848,  646,  464  nm 


The  second  sample  case  presents  a  subset  of  a  SPOT-5  scene  from  Rumania,  acquired  22  May  2005, 
containing  a  large  percentage  of  cloud  areas  (Fig.  4).  A  lot  of  details  can  be  seen  in  the  de-shadowed 
image  that  are  hidden  in  the  original  scene. 
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Fig.  4  SPOT-5  scene,  RGB=836,  657,  555  nm.  Courtesy  European  Space  Imaging. 


The  third  example  shows  de-shadowing  of  building  shadow  areas  from  a  Munich  scene  recorded  by 
Ikonos,  17  Sept.  2003,  see  Fig.  5.  Although  the  contrast  in  the  de-shadowed  image  is  drastically  reduced 
as  it  should  be,  the  zoomed  view  in  shadow  areas  reveals  additional  information  that  cannot  be  found  in 
the  original  scene. 


Fig.  5  Ikonos  scene  of  Munich  center,  RGB=802,  665,  551  nm. 

Top  left:  original,  right:  after  de-shadowing.  Courtesy  European  Space  Imaging. 
Bottom  shows  a  zoomed  view:  original,  shadow  map,  de-shadowed  image. 


4.  SUMMARY 

A  fully  automatic  de-shadowing  method  for  spectral  imagery  has  been  developed.  An  essential  part  of  the 
technique  is  the  concept  of  core  shadow  and  expanded  shadow  regions  to  reduce  the  number  of  shadow- 
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pixel  misclassifications  significantly.  The  algorithm  was  implemented  with  three  tunable  parameters:  the 
size  of  the  core  shadow  mask  (small,  medium,  large),  the  maximum  range  of  the  unsealed  shadow 
function,  and  the  minimum  fractional  direct  solar  illumination  of  the  darkest  shadow  areas,  the  depth  of 
shadow.  Once  the  parameters  have  been  set  the  algorithm  runs  fully  automatic.  Default  values  are 
provided  that  usually  yield  good  de-shadowing  results.  However,  for  optimum  performance  these 
parameters  might  have  to  be  tuned  depending  on  the  scene. 
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